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Branching of bilayer membranes appear in the inverted hexagonal phase as weU as in metastable 
states of the lamellar phase such as membrane fusion intermediates. A method for estimating the 
line tension of the branching junction is proposed for molecular simulations. The line tension is 
calculated from the pressure tensor of equiangular ly branched membranes. The simulation results 
agree very well with the theoretical prediction of Hamm and Kozlov's tilt model. The transition 
between the lamellar and inverted hexagonal phases is also investigated using the tilt model. 



I. INTRODUCTION 

Amphiphilic molecules, such as lipids and detergents, 
self-assemble into various structures: spherical or cylin- 
drical micelles- lamella; and hexagonal, cubic, and sponge 
structures [mtM- The bilayer membrane in the lamellar 
phase is the basic structure of the plasma membrane and 
the intracellular compartments of living cells. It is known 
that biological membranes contain a substantial amount 
of lipids that can adopt the inverted hexagonal Hn phase 
upon isolation 0-0] • Thus, it has been speculated that 
non-bilayer structures play some roles in biological func- 
tions such as membrane fusion and other membrane con- 
tact phenomena such as tight junctions. 

Several mechanisms of membrane fusion have been pro- 
posed Among these, the stalk model is widely 
accepted. According to this model, two vesicles form a 
metastable structure, the stalk intermediate, where the 
outer monolayers (leaflets) are connected with a cylin- 
drical (stalk) structure. Next, radial expansion of the 
stal k p , [l3| or a pore opening on the side of the stalk 
O leads to either complete fusion or the formation of 
a hemifusion diaphragm (HD) intermediate. In the HD, 
the inner monolayers of two vesicles form a bilayer mem- 
brane, which connects with two bilayer membranes of 
the vesicles via a branching junction. Many experiments 
support that fusion proceeds through a metastable inter- 
mediate state, which is considered to be the stalk or HD 
intermediate. Recently, molecular simulations of mem- 
brane fusion and fission have been reported by several 
groups [TTI. [T2|. IT^ - [23| . In most of these simulations, the 
stalk intermediate is formed first. However, the path- 
way of the stalk to complete fusion is not unique and 
is dependent on the simulation conditions. Thus far, it 
has not been determined which membrane properties and 
external conditions determine the fusion pathways. It is 
important to estimate the free energy of each fusion in- 
termediate to understand the fusion pathways. 

Since the transition between lamellar and Hu phases 
requires the same type of topological changes as mem- 
brane fusion, it is considered to proceed through simi- 
lar pathways [H, HI] . Recently, the phase of a hexago- 
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nally ordered stalk structure has been observed between 
lamellar and Hjj phases [25|. Here, we consider the 
Hii structure with a long cell length. The bilayer mem- 
branes are connected with hexagonally arranged branch- 
ing junctions. When the branching junctions of bilayers 
are isolated, i.e., sufficiently far from other junctions, the 
free energy of a junction is proportional to the junction 
length. Thus, it can be estimated as a line tension. In 
this paper, we propose a method for estimating this line 
tension for molecular simulations. This tension quantifies 
the stability of the HD fusion intermediate. 



To simulate lipid membranes on long and large scales, 
various coarse-grained molecular models have been pro- 
posed (see review articles (26l - l29| ). Recently, the bottom- 
up approaches have been intensively investigated to con- 
struct coarse-grained molecular models, in which poten- 
tial parameters are tuned from atomistic simulations [1- 
1131 • We, on the other hand, choose an opposite top-down 
approach to construct the coarse-grained molecular mod- 
els, in which potentials are based on the continuum the- 
ory. Here, we employ a spin lipid molecular model [SSf, 
which is one of the solvent-free molecular lipid models 
j28| . In this model, the membrane properties (bending 
rigidity, line tension of membrane edge, area compression 
modulus, lateral diffusion coefficient, and flip-flop rate) 
can be varied over wide ranges. 



In Sec. ini the simulation method and the membrane 
properties of the simulation model are described. In Sees. 
Illll and irvl the line tension is estimated by the simulation 
and the tilt model, respectively. The line tension is cal- 
culated from the pressure tensor of hexagonally arranged 
membrane branches. Previously, the line tension of the 
membrane edge has been estimated from the pressure 
tensor of a membrane strip [33 - l36| . For branching junc- 
tions the membrane can hold finite surface tension, unlike 
the strip, hence, the membrane area should be carefully 
treated. The tilt model proposed by Hamm and Kozlov 
[37I [HI reproduces the simulation results very well. In 
Sec. |Vl the rupture dynamics of the branches and the 
HD intermediates of the vesicles are shown. Finally, a 
summary is given in Sec. IVIl 
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II. SIMULATION MODEL AND METHOD 
A. Membrane Model 

A spin lipid molecular model [35| is used to simulate 
bilayer membranes. Each (i-th) molecule has a spherical 
particle with an orientation vector u^, which represents 
the direction from the hydrophobic to the hydrophilic 
part (|ui| = 1). There are two points of interaction in 
the molecule: the center of a sphere r| and a hydrophilic 
point r° = rf + Uia. The molecules interact with each 
other via the potential, 
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where = rj 



5 1 — ^ij — ^ij/^ij^ S-nd k^T 

is the thermal energy. The molecules have an excluded 
volume with a diameter a via the repulsive potential, 
Urepir) = exp[— 20(r/CT — 1)], with a cutoff at r = 2Aa. 

The second term in Eq. (HI) represents the attractive 
interaction between the molecules. An attractive multi- 
body potential Uatt{pi) is employed to mimic the "hy- 
drophobic" interaction. This potential allows the forma- 
tion of the fluid membrane over wide parameter ranges. 
Similar potentials have been applied in other membrane 
models [39l - l4l| and a coarse-grained protein model [i^ . 
The potential C/att(Pi) is given by 

U^ttifr) = 0.251n[l +exp{-^4(p, - p*)}] - C, (2) 

with C = 0.25 ln{l -I- exp(4/9*)}. The local particle den- 
sity Pi is approximately the number of particles r| in the 
sphere with radius ratt- 



where /cutC?") is a C°° cutoff function, 
fcnt{r) 



exp{A(l+ (r<reut) 

(r > Teut) 



(3) 



(4) 



with n = 6, A = ln(2){(reut/ratt)" - 1}, ^att = 1.9a 
(/cut('"att) = 0.5), and the cutoff radius rcut — 2.4cr. The 
density p* in f7att(Pi) is the characteristic density. For 
Pi < p* — l, Us^ttifi) acts as a pairwise attractive potential 
while it approaches a constant value for pi > p* + 1. 

The third and fourth terms in Eq. ([1} are discretized 
versions of the tilt and bending potentials of the tilt 
model j37l.[38j, respectively. A smoothly truncated Gaus- 
sian function is employed as the weight function 
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with n = 4, Tga = 1.5(7, and Tcc = 3(t. All orders of 
derivatives of /cut('') and Wyaisi"!^) are continuous at the 
cutoff radii. The weight is a function of to avoid the 
interaction between the molecules in the opposite mono- 
layers of the bilayer. If rf^- is employed instead, a single- 
layer membrane is formed [43l |. 



B. Simulation Method 

The bilayer membrane is simulated in the NVT en- 
semble (constant number of molecules N , volume V , and 
temperature T) and Brownian dynamics (molecular dy- 
namics with Langevin thermostat) is employed. The mo- 
tion of the center of the mass rp = (r^ -I- r°)/2 and the 
orientation are given by underdaniped Langevin equa- 
tions: 



dt 



m- 
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dt 

dt 

= -C,u;,-H(gKi)+f, 
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where m and / are the mass and the moment of iner- 
tia of the molecule, respectively. The forces are given 
by fp = -dU/drf and f[ = -dU/du, with the per- 
pendicular component = a — (a • Ui)ui and a La- 
grange multiplier Al to keep uf — 1. According to 
the fluctuation-dissipation theorem, the friction coeffi- 
cients Cg and Cr and the Gaussian white noises gp(t) 
and gj(t) obey the following relations: the average 
(gtA^)) = and the variance (fff (^Osf., (^2)) = 
S{ti - t2), where ai, 02 G {x, y, z} 
and /3i,/32 G {G,r}. The Langevin equations are inte- 
grated by the leapfrog algorithm [H, Hj] ■ 

The results are displayed with a length unit of a, an 
energy unit of /cbT, and a time unit of tq — Qqu^ / k^T . 
m = CqTq, I = CrTo, Cr = CGa^ Ai = 0.005ro, and p* = 
14 are used. The error bars of the data are estimated 
from the standard deviations of three or six independent 



C. Membrane Properties 

The potential-parameter dependence of the membrane 
properties (bending rigidity k, line tension F of mem- 
brane edge, area compression modulus, lateral diffusion 
coefficient, and flip-flop rate) are investigated in detail in 
our previous paper ,35]. In this model, these properties 
can be varied over wide ranges. The line tension F of 
the membrane edge can be controlled by varying e and 
Cbd- The membrane has a wide range of fluid phases, 
and the fluid-gel transition point can be controlled by 
p* . The area compression modulus Kp^ can be varied by 
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The flip-flop rate can be varied by kti\t and k 
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FIG. 1: Bending rigidity k dependence on (a) fctiit and (b) 
fctiit + fcbond for e = 2 and Cbd = 0. The lines with squares, 
circles, and triangles represent data for fcbend = 0, fcbcnd = 
fctiit, and fcbond = 8, respectively. The crosses and diamonds 
represent data for fctiit = 4 and 8, respectively, for several 
values of fcbend- The dashed line shows fv/fer = 2.2(A:tiit + 

fcbcnd) + 6. 

The bending rigidity k is linearly dependent on fctut 
and fcbend- The spectra of undulation modes = 
k-QT/Kq^ are widely used to estimate k, for tensionless pla- 
nar membranes. In Ref. [s^l, we use a least-squares fit to 
1/ {\h{q)\'^) — nq^ /k-QT for low frequency modes q < (/cut 
with ((7cut/7i')^ = 0.015. More recently, we have found 
that the extrapolation to ^cut yields a better k esti- 
mation [4^. Re-estimated values of k are shown in Fig. 
[T] for the tensionless membrane using this extrapolation 
method. Slig htly higher values are obtained and the slope 
2 in Ref. [33 becomes 2.2 as re/fceT — 2.2(fctiit + fcbond)+6 
for e = 2 [see Fig. ^h)]. 

The bending elasticity generated by the bending and 
tilt potentials can be derived from the continuum theory 
[i^Ii^ as discussed in Refs. [H,!!^. When the orienta- 
tion vectors are equal to the normal vectors n of the 
membrane without tilt deformation, the bending and tilt 
energies of the monolayer membrane are given by 

= JdA ^Scnd + <ilt (g^ ^ - Cof 

-(«bc„d + <ilt)ClC2 (9) 

in the continuum limit, where Ci and C2 are the two prin- 
cipal curvatures of the monolayer membrane. The bend- 
ing rigidity and the saddle-splay modulus of the mono- 
layer are given by Kmo = ^bond + '^tiit and = -k, re- 
spectively. The spontaneous curvature is given by Cq = 




FIG. 2: Snapshot of the bilayer membranes with hexagonally 
arranged branches for A'^ — 3072, Lx = VSLz/2 = 5Qa, Ly = 
16.94(7, e = 2, fcbend = 8, fctiit — 8, and Cbd ~ 0. (a) Periodic 
images, (b) Magnified snapshot of the branching junction. 
The red and yellow himispheres represent the hydrophilic and 
hydrophobic parts of the molecules, respectively. 



{'«bond/('*bond + '«t*iit)}C'bd/''nb, whcrc r„b is the nearest- 
neighbor distance. By assuming a hexagonal packing of 
the molecules, the bending rigidities generated by the 
bending and tilt potentials are estimated as '^bcnd/^B-^ ~ 
V3fc bcnd'^cv (r„b) and K^jj^/fceT = y3fctiitWcv(?'nb)/2, re- 
spectively. Thus, the bending rigidity k, of the bilayer 
membrane is estimated as k = 2Kmo — (2.1fcbcnd + 
l.lfctiit)fcB2^ from fnb — 1.05(7 and u'cv(l-05(7) — 0.61. 
This estimation of k shows quantitative agreement with 
the simulation results for fcbond, while the prefactor of fctiit 
is half that of the numerical estimation. Deviation from a 
hexagonal lattice and tilt deformation likely change the 
prefactor of fctiit. Using — 2k^j[^, the spontaneous 
curvature Cq ~ {fcbond/ (fcbond + fctiit ^jCbd/o- is obtained, 
which agrees well with Co calculated from a membrane 
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FIG. 3: Dependence of the junction line tension Fbr and sur- 
face tension 7 on the membrane area a^y per molecule for 
£ = 8, fcbcnd = 8, fctiit = 8, and Cbd = 0. (a) The line ten- 
sion Fbr calculated from Eq. (|12p with A^b ~ axyN/2 (o, x) 
and j4mb = '2xLy (□,<>). The symbols (x,o) and (o,n) rep- 
resent data for A*' = 1536 and = 3072, respectively, (b) 
The molecular area a^y for the bilayer in the middle of the 
simulation box is compared with the geometrically estimated 
area ol = AL^Ly/N for N = 1536 (x) and iV = 3072 (o). 
(c) The surface tension 7 calculated from the branched mem- 
branes (o) for A*' = 3072 and the previous simulation of the 
planar bilayer membrane (A) for A'^ = 512. (d) The differ- 
ences in the values of 7 estimated from the three methods for 
A'^ = 3072. The circles (o) and squares (□) represent 7^ — "jx 
and 7mid — 7a:, respectively. The tensions "fx and 72 are cal- 
culated from the pressure tensors P^x and Pzz, respectively. 
The tension 7inid is calculated from the middle bilayer. 



tube and strip [43|. Positive spontaneous curvature in- 
dicates that the hydrophilic head is larger than the hy- 
drophobic tail of amphiphilic molecules. 



III. MEASUREMENT OF LINE TENSION OF 
BRANCHING JUNCTIONS 

To measure the stress of branched membranes, we con- 
sider hexagonally arranged branches (see Fig. [5]). These 
have an inverted hexagonal structure with a long cell 
length. Three membranes equiangularly connect with 
each other at the branching junction. The periodic 
boundary along the x direction is shifted as shown in Fig. 
[21a) . Normal periodic boundary conditions are employed 
in the y and z directions. Thus, the periodic images are 
given by {xi + n^Lx, yi + riyLy, Zi + n^Lz + n^L^JI) = 
{xi,yi,Zi), where Ux, riy, and Uz are arbitrary inte- 
gers for a simulation box with side lengths L^, Ly, and 
Lz. Length Ly is varied to change the membrane area 
and surface tension. The other side lengths are fixed 
as Lx = V3Lz/2 = 28a and 56a for N = 1536 and 
3072, respectively. When the membrane thickness is ne- 
glected, the total membrane area in the simulation box is 
— 2LxLy, since the area of one membrane is 2LxLy/S 
and three membranes exist in the box. The length of one 
branching junction is Ly, and two junctions are in the 
box. The origin {x,y,z) = (0,0,0) is set at the center 
of the simulation box and is kept at the center of the 
horizontal membrane. 

The line tension of the branching junction can be cal- 
culated from the diagonal components of the pressure 
tensor 

dU 

Pea = {NkBT - a^)/^' (10) 

where a G {x,y,z} and V — LxLyL^- The stresses 
should be generated by the line tension Fbr of the branch- 
ing junction, the surface tension 7, and the bulk pressure 

^bulk- 

P.xV = PbuikV - ^LxLy, (11) 

PyyV = PbulkV^-7Anb-2FtaXy, (12) 

P.zV = Pbuikl^-^7iyi.- (13) 

In the membranes with an angle of ±7r/3 to the xy plane, 
the surface tension yields stresses 7/2 and (•\/3/2)7 in the 
X and z directions, respectively. Since the membranes 
can be moved in the x and z directions, these stresses 
are balanced i.e. Pxx = Pzz- Indeed, Eqs. pT|) and 
give the same value. Since all membranes are parallel to 
the y direction, the total membrane area ^mb in the sim- 
ulation box is used in Eq. p^ . Since the critical micelle 
concentration (CMC) is low and no isolated molecules 
appear, the bulk pressure Pbuik is negligibly small in our 
solvent-free simulations. 

To estimate the line tension Fbr, the surface tension 
7 has to be calculated first. In the simulations, 7 can 
be estimated from Pxx or P^z- We called them ^x and 
7^, respectively. We also calculated 7 from the mid- 
dle of the bilayer membrane which is parallel to the 
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FIG. 4: Dependence of the line tension Fbr of the branching 
junction on (a) ktut and (b) Cbd- (a) The squares, circles, 
and triangles represent data for fcbend ~ 0, fcbcnd ~ fctiit, and 
fcbcnd = 8, respectively, at e = 2 and Cbd = 0. (b) The 
symbols (o, □) and (x,o) represent data for e = 2 and 8, 
respectively, at fcbcnd ~ fctiit = 4 or 8. The solid lines in (a) 
and (b) are obtained from Eq. (|22p . 



xy plane using the definition for the planar interface, 
{Pxx+Pvv)/'2)L, for -L^/12 <x< L^/U. 
The estimated values of 7 from these three methods agree 
very well [see Fig. [31Jd)]. We use the average value 
7 = ilx + 72)72 as the surface tension to reduce sta- 
tistical errors. Here, we only consider membranes with 
small or positive surface tensions, i.e., exclude the buck- 
led membranes under lateral compression, since the sur- 
face tension becomes anisotropic in buckled membranes 

Next, we consider the membrane area. The pro- 
jected area a^y per molecule is calculated from the mid- 
dle of the membrane at —Lx/12 < x < Lx/12 as 
Gxy = (La;Lj;/6)/(iVmid/2) where Nyaid is the number 
of the molecules at —Lx/12 < x < Lx/12. The a^y 
dependence of 7 calculated from branched membranes 
agrees with that obtained from planar membranes [see 
Fig. [Stc)]. This slightly smaller than the area 

ol = 2LxLy/{N/2) calculated from the geometry. The 
decreases are 4% and 2% for TV — 1536 and 3072, re- 
spectively, as shown in Fig. [31[b). These deviations are 
caused by the membrane thickness. In Fig. [Sj the length 
of the membrane surface is shorter than the center line 
of the membrane. Thus, the monolayer (leaflet) area of 
the bilayer membrane becomes smaller than ol. 



is estimated from the pressure tensor and a^y, 

T^.^{{Pxx+P..)^~Pyy]^. (14) 



In explicit solvent simulations (Pbuik 7^ 0), 7 or Pbuik 
should be calculated from a local region in the simulation 
box. When ^i-aid is used, the line tension is given by 



Fbr - (— 



p I V' 7mid-^2; 2axy 

^^I2l7y 2 



.( \— ^""^f ^ 



(15) 

for explicit solvent systems. When Pbuik is estimated 
locally, the line tension is given by 



Fbr = { (Px 



P. 



s ) Pyy ^~ -^bulk 
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The tension Fbr is almost independent of 7 [see Fig. 
inija)] . The estimated values of Fbr for two system sizes 
N = 1536 and 3072 agree very well, so that neighbor- 
ing junctions are sufficiently far away to avoid the fi- 
nite size effects of the membranes. If we naively use 
ylmb = 2LxLy, the obtained values of Fbr decrease with 
increasing 7, and even become negative [see Fig. [Sja)]. 
Thus, although the area difference between axy and ol is 
very small, it has a large influence on the estimation of 
Fbr. The monolayer membrane area a^yf^ 12 should be 
employed. 

We investigated the Fbr dependence on the potential 
parameters in our membrane model, as shown in Fig. 
The tension Fbr is roughly linearly dependent on /ctiit and 
Cbd while it is almost independent of e. These dependen- 
cies can be understood by the energy minimization in the 
tilt model as explained in the next section. 



IV. TILT MODEL 

We analyze the line tension of the branching junction 
and the stability of the inverted hexagonal Hn phase 
using the tilt model proposed by Hamm and Kozlov [13, 
[38j . The orientation vector u of the lipid molecules can 
deviate from the normal vector n of the dividing surface 
of the bilayer. The tilt vector is defined as t = u/u-n— n, 
which is parallel to the dividing surface (see Fig. [S]). The 
tilt tensor is the derivative of the vector t: = dta/dp, 
where a, /3 € {x,y} for a membrane parallel to the xy 
plane. In the tilt model, the free energy of monolayers of 
the tensionless membranes is given by 



F 



2 ^ " 



Co)' 



3det(i2) + 



(17) 

where and det(t|) are the trace and determinant of 
the tilt tensor, respectively. The coefficients Kmo and 
Finally, the line tension Fbr of the branching junction are the bending rigidity and saddle-splay modulus of the 
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FIG. 5: Schematic representation of tiie tilt model. 
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FIG. 6: The cell length L (a) and free energy F (b) of inverted 
hexagonal Hu phase with respect to the lamellar phase. The 
solid and dashed lines are obtained from the tilt model with 
variation calculus [Eqs. (|21[l ] and the linear approximation 
[Eq. (O], respectively. 



monola yer , respectively, in the Helfrich-Canham bending 
energy |45l. l46j . The coefficient Kg is the tilt modulus. 

First, we consider the free energy of a horizontal mono- 
layer in Fig. [5] Here, we use the same coordinates as the 
simulations [n — (0,0,1)]. Since the molecules tilt only 



in the x direction. 



ty = 0, the free energy of the 



monolayer is expressed as 



hex 



2L„ 



L/2 



dx, 
(18) 

The boundary conditions are ^^^(O) ~ 
tan(7r/6) — — 1/\/3, where L is the 
side len gth of the hexagonal cell, L = 2Lx/i. Hamm and 
Kozlov ^31 obtained 



where t'^ = dtx/dx 
and t^{L/2) = - 



Fhc 



A 



Co 



18' 



(19) 



from the assumption that t'^ is constant. Here, we use 
tx{x) which minimizes -Fhcx, instead of the constant t'^ 
assumption. The calculus of variation gives t'^{x) = 



{Kg/Kyno)tx{x). Hence, the molecules are tilted by 

sinh(a::/A) 



txix) 



V3sinh(L/2A)' 



(20) 



where A = y^^moT'^ is the characteristic length. Then 
the free energy is obtained as 



hex 



sinh(i/A) 



imoLy 3A(cosh(i/A) 
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2Co 
V3 



Co^L 



(21) 



When the distance L between the junctions is sufficiently 
large {i.e. L ^ A), the first term on the R.H.S of Eq. 
(PT|) becomes 1/3A. 

The free energy Fiam of the lamellar phase is given 
by Fiam/^ = K,„oCo/2- The line tension of the branch- 
ing junction is the energy required to form a junction 
per unit length, Tbr^y = 3(Fhcx — -Flam), since one junc- 
tion links six monolayers and one monolayer connects two 
junctions. Thus, Tbr is given by 



Fbr = y/KmoKeil + 2V3ACo), 



(22) 



for L > A. 

Next, we compare the line tension of the tilt model 
with our simulation results. Molecular tilts are seen 
around the junctions in simulation snapshots [see Fig. 
[5Jb)]. In our membrane model, the bending rigidity 
and spontaneous curvature of the monolayer are given by 

Kmo/k^T ~ l.l(fcbond + ^tilt) and Co ^ {Kbcnd / (ti^b cnd + 

Ktiit)}C'bd/f, respectively, as described in Sec. Ill CI The 
tilt modulus Kg is estimated as follows. When thermal 
fluctuations are neglected, the tilt energy Utat of a flat 
membrane with a constant tilt is given by 



-Ftilt ^tilt, 2 AT AT / \ 

kbT 4 



(23) 



where N-a^ is the mean number of neighboring molecules. 
Here the average ((t-f)^) = t^/d is used for isotropic unit 
vectors f with a spatial dimension d = 2. In the contin- 
uum limit, i^tiit should correspond to -Fhox = ugAtx'^ /2 of 
the tilt model. Thus, Kg is given by Kg ~ 2ktiitkBT / Uxy 
from A^'nb = 6 and (wcv) = u;cv(l-05f7) = 0.61. The line 
tension Fbr calculated from Eq. reproduces the sim- 
ulation results very well for most of the parameter ranges 
in Fig. m The tilt model overestimates Fbr for fcbcnd = 
in Fig. Ufa). This deviation is likely caused by a too 
small characteristic length A. It is less than the molec- 
ular diameter cr: A ~ 0.8a and l.la at fcbcnd = and 
fcbcnd = fctilt, respectively, for the tensionless membranes 



with 



1.2a2. 



Finally, we discuss the stability of the inverted hexago- 
nal Hjj phase using the tilt model. For — ACq < 1/2 \/3, 
the free energy per area Jhcx/^ monotonically decreases 
with increasing cell length L [see Eq. (I?T|) ]. In this re- 
gion, the lamellar phase {L — )■ oo) is more stable than the 
Hjj phase. For —ACq > l/2-\/3, -Fhcx/^ is minimum at 
flnite L and the Hjj phase is more stable than the lamel- 
lar phase. In this region, the line tension Fbr is negative. 
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FIG. 7: Sequential snapshots of the membrane rupture for 
iV = 3072, £ = 2, fcbend = 8, fctiit = 8, and Cbd = 0. (a) 
t/ro = 2555. (b) t/ro = 2565. (c) I/tq = 2570. (d) t/ro = 
2575. (e), (f) t/ro = 2580. 



and the Hu hexagonal structure grows. Figure[6]shows L 
and (^hex — -Fiam)/^ obtained by the present method and 
the previous linear approximation 37|. Under the linear 
approximation, (Jhex--Fiam)M = '«e/18-K„ioC'o/2 with 
L = — 2/\/3Co, so that the Hu phase has a lower energy 
for — ACo > 1/3. As — ACq increases, the two lines in Fig. 
[6] approach each other. Thus, the linear approximation 
works well in the Hjj phase at — ACq > 0.5, while not for 
-ACo < 1/2 V3. 



V. DYNAMICS 



A. Rupture 



1 1 1 1 1 1 1 1 1 1 r 
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FIG. 8: Parameter e dependence of (a) surface tension 7cut of 
the membrane rupture and (b) the line tension of the mem- 
brane edge for fcbend = fctiit = 8. The triangles, circles, and 
squares represent data for Cbd = —0.2, 0, and 0.2, respec- 
tively. 



point. As £ increases, 7rup increases linearly as shown 
in Fig. 13a), although Fbr is independent of e. The line 
tension F of the membrane edge also increases but it is 
not linear with respect to e [see Fig. Hl^b)]. 

After the membrane rupture, a branched membrane 
becomes a line of a membrane edge and a bent membrane. 
The free energy of the former and latter states of the 
tensionless membranes are given by F^^/Ly = Fbr and 
Fiup/Ly = F -I- (7r/6)K/i?, respectively, where R ^ 5a is 
the curvature radius of the bent membrane. Thus, the 
branched membrane is kinetically stable for the rupture 
when Fbr < T for 7 = 0. Furthermore, even when i^br > 
i^rup) the branched membrane can remain in a metastable 
state to overcome the free energy barrier to rupture. In 
the present spin molecular model, Fbr and F can be varied 
separately by ktm or /cbend and e, respectively. 



A membrane rupture occurs under sufficiently high 
surface tension 7. We investigated the rupture of 
branched membranes as Ly is increased slowly at a rate 
dLy/dt = 0.0002(t/t. We checked that this rate is suffi- 
ciently small. No difference is detected when a rate twice 
as fast, i.e., dLy/dt — 0.0004o-/r, is used. Figure |7| shows 
the rupture of the branched membranes. A pore opens 
on the side of the branching junction and then expands 
along the junction. Since the junction is less stable than 
the middle of the bilayer membranes, the rupture always 
occurs on the side of the branching junction. 

The surface tension 7rup at the rupture is estimated 
from the extrapolation of the Ly-j line to the rupture 



B. Fusion Intermediate 

The branching junctions appear in the hemifusion di- 
aphragm (HD) intermediate in membrane fusion [see Fig. 
^a)] . The HD intermediate was proposed in the original 
stalk model of membrane fusion '48] and was observed in 
experiments [iQl and in some of the molecular simulations 
[lE[li,|2l|,|2j7^Here we focus on the structure of the HD 
intermediate and the detail of the fusion dynamics of the 
present model will be described somewhere else. Follow- 
ing the procedure in our previous fusion simulations jl7| , 
a spherical particle with a radius of 4(t is placed at the 
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FIG. 9: Branched structures in partially fused vesicles for 
£ = 8, fcbend ~ 8, fctilt ~ 8, and Cbd = 0. Sliced snapshots 
are shown in (a) two and (b) three hemifused vesicles for 
= 4000 and A'^ = 6000, respectively, (c) All molecules are 
shown for the snapshot in (b). 

center of each vesicle and constant external forces are ap- 
plied to the spherical particles. The membranes pinched 
by the particles are fused either directly or via the HD 
intermediate. The fusion dynamics are qualitatively sim- 
ilar to our previous simulations [l^. When the HD is 



formed, we removed the spherical particles and then re- 
laxed the vesicles to a stable structure. The resultant 
structures are shown in Fig. [9l 

The HD of two vesicles has a circular branching junc- 
tion. As the radius i?br of the branching junction in- 
creases, the bending energy is reduced because the flat 
membrane area also increases. However, the energy of 
the junction increases as 27ri?brrbr. Thus, the radius 
is determined by the balance of these energies. Figures 
ini (b) and (c) show the HD of three hemifused vesicles. 
Four lines of the branching junction meet at the centers 
of the front and back membranes. These vesicle shapes 
resemble the shape of soap bubbles. In the bubbles, com- 
petition between surface tension and pressure maintains 
their shapes, while in the vesicles, the bending energy 
and line tension play these roles. 



VI. SUMMARY 

We have studied the branching of bilayer membranes. 
The free energy of the branching junction can be treated 
as line tension when the distance between the junctions is 
sufficiently large. The line tension Fbr of the junction is 
calculated from the pressure tensor in solvent-free molec- 
ular simulations. The simulation results agree well with 
the prediction of the tilt model. 

The stability of the branching membranes depends on 
the boundary conditions. High surface tension induces 
membrane rupture on the side of the branching junction. 
In the hcmifusion diaphragm intermediate, the branch- 
ing junctions can be formed by the balance of the bend- 
ing energy and line tension. Several kinds of proteins 
are known to modify membrane curvatures in living cells 
[i^ - fsij . The branched structures induced by the jibsorp- 
tion of a colloid [15|], DNA gli], and peptide ^] have 
been reported in molecular simulations. The branching 
structures of the biomembranes may be stabilized or con- 
trolled by the absorption or insertion of proteins to the 
membranes in living cells. 

The self-assembled structures of amphiphilic molecules 
can be clarified from the packing parameter fH of am- 
phiphilic molecules or the spontaneous curvature of the 
monolayer. However, these quantities are not easy to 
measure. Alternatively, the line tension of the branching 
junction may be suitable to quantify the tendency of the 
amphiphilic molecules to form inverted structures. 
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